Flame Enhancement and Quenching in Fluid Flows 
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We perform direct numerical simulations (DNS) of an advected scalar field which diffuses and 
reacts according to a nonlinear reaction law. The objective is to study how the bulk burning rate 
of the reaction is affected by an imposed flow. In particular, we are interested in comparing the 
numerical results with recently predicted analytical upper and lower bounds. We focus on reaction 
enhancement and quenching phenomena for two classes of imposed model flows with different ge- 
ometries: periodic shear flow and cellular flow. We are primarily interested in the fast advection 
regime. We find that the bulk burning rate v in a shear flow satisfies v ~ aU + b where U is the 
typical flow velocity and a is a constant depending on the relationship between the oscillation length 
scale of the flow and laminar front thickness. For cellular flow, we obtain v ~ U 1 ' 4 . We also study 
flame extinction (quenching) for an ignition-type reaction law and compactly supported initial data 
for the scalar field. We find that in a shear flow the flame of the size W can be typically quenched 
by a flow with amplitude U ~ aW. The constant a depends on the geometry of the flow and tends 
to infinity if the flow profile has a plateau larger than a critical size. In a cellular flow, we find that 
the advection strength required for quenching is U ~ W 4 if the cell size is smaller than a critical 
value. 

PACS numbers: PACS numbers: 47.70.Fw, 47.27.Te, 82.40.Py 



I. INTRODUCTION 

Turbulent combustion in premixed flows is a widely 
studied topic in both scientific and industrial settings (see 
e.g. Q |2|). The interest in the subject is due to 
an important influence that advection can have on the 
reaction process: both experimental [O, |4^| and theoret- 
ical §, g, H, H, 11 H, S ill work shows that 
the propagation speed of the flame can be significantly 
altered by the fluid flow. Specifically, moderately intense 
levels of turbulence have the tendency to accelerate the 
flame speed v beyond its laminar value v . The mecha- 
nism and the extent of the flame acceleration depend on 
the particular regime of burning The general rea- 

son for the enhancement is that the fluid motion distorts 
the flame front, increasing the reaction area. On the 
other hand, if the advection is too strong, it can lead to 
the flame extinction. The critical strength of advection 
which leads to quenching depends on the extent of the 
flame, strength of reaction and diffusion, and properties 
of the flow. 

At this stage, it is unreasonable to expect a complete 
analytical theory describing the process of combustion in 
a fluid phase. Indeed, detailed modelling of the phenom- 
ena involves solving a reaction-diffusion system involving 
temperature (or energy) and concentrations of reactants 
coupled with compressible Navier-Stokes equations de- 
scribing motion of the mixture p9| , j57f . Therefore, most 
of the studies in this field which seek analytical conclu- 
sions use heuristic reasoning or simplified models, which 
may approximately describe the system in certain com- 
bustion regimes. Some of the combustion regimes are 
relatively well-understood, such as the so-called flamelet 
regime, where flame thickness is small compared to the 



fluid velocity scales. The geometric optics approximation 
where the propagation of the front is ruled by Huygens 
principle is often used as a starting point in the analysis 
of this regime (see e.g. |34|, |38|). 

Our goal here is to study one of the most widely used 
PDE models of combustion, namely the scalar reaction- 
diffusion equation with passive advection: 
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Here T is the normalized temperature, < T < 1, u is 
the fluid velocity, which we assume is incompressible, n 
is thermal diffusivity, and r is the typical reaction time. 
In the absence of fluid velocity Eq. (|l|) admits flat prop- 
agation front with laminar burning velocity of the order 
of v ~ \/k/t and characteristic thickness of the order 
of S ~ s/kt. The model (0) can be derived from a more 
complete system under assumptions of constant density 
and unity Lewis number (the ratio of material and tem- 
perature diffusivity), as shown, for instance, in |Tfl |. The 
equation ([!]) has a more general applicability than the 
geometrical optics approximation; moreover, as we will 
discuss below, the geometrical optics limit can be ob- 
tained from (|l|) in a certain parameter range. 

We will consider reaction rates R(T) of two types, KPP 
(Kolmogorov, Petrovskii, Piskunov) [^7j and igni- 
tion. The KPP type is characterized by the condition 
that the function R(T) is positive and convex on the in- 
terval < T < 1. This reaction type is used often in 
problems on population dynamics (see e.g. [^, |26|), but 
is relevant in combustion modelling, for example in some 
autocatalyctic reactions |28| . A reaction term of ignition 
type is characterized by the presence of critical ignition 
temperature, such that the function R(T) is identically 
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zero below ignition temperature. This type of reaction 
term is used widely to model combustion processes (see 
e.g. jl9|, H^]), in particular approximating the behavior of 
Arrhenius-type chemical reactions which vanish rapidly 
as temperature approaches zero. 

Our main goal is to gain insight into the question of 
how the geometry and the amplitude of the fluid flow 
influence the combustion process. Our study is partly 
motivated by recent analytical work j(| 19, 



where rigorous bounds on combustion enhancement and 
quenching are proved. We test the sharpness of results in 
H |l^, |2^, [35| , and in addition derive new predictions. 
We consider two classes of flows. The first is shear flows, 
a representative of a wider class of flows, called "perco- 
lating" in Jl9| , which have open streamlines connecting 
distant regions of the fluid. The second class is cellular 
flow, where the streamlines are closed and the flow con- 
sists of isolated cells. For each class of flows, we study 
both flame enhancement and quenching. 

For flame enhancement study, we consider initial tem- 
perature in the form of the laminar front, with T = 1 in 
the semi-infinite region behind the front and T — in the 
the semi-infinite region ahead of the front. Distorted by 
imposed flow, the flame front propagates as a travelling 
wave with velocity higher than laminar. The goal of the 
flame enhancement study is to obtain relations between 
flame propagation speed v and the properties of the flow, 
especially for the large advection velocities. 

In the case of quenching phenomena, we consider ini- 
tial temperature to be non-zero in a finite region. Since 
quenching cannot occur for the KPP-type source term 
p3| , ffij , we use the ignition- type reaction term. As shown 
by Kanel, there is a critical size W Q of initial hot re- 
gion below which the flame will be extinguished by diffu- 
sion alone, e.g. with no advection, when the temperature 
drops below the threshold and reaction ceases before the 
flame establishes a steady travelling wave configuration 
p2| . When advection is present, the fluid flow stretches 
the initially hot region so that it can be quenched by 
diffusion; hot regions of the size much larger than W a 
can be quenched in this manner. Our goal has been to 
understand how the geometry and amplitude of the flow 
influence the size of the band W of the initial hot region 
that can be quenched. 



II. NUMERICAL SETUP AND METHOD 

The simulation is set in two space dimensions, in a ver- 
tical strip of width L with periodic boundary conditions 
in x direction (Fig. [j]). In reaction enhancement studies, 
the initial temperature was set to T = 1 in the lower half 
of the domain and to T = in the upper half of the do- 
main. In the quenching studies the initial temperature 
was set to T = 1 in a horizontal band of width W in 
the center of the domain, and to T = elsewhere. The 
interfaces between hot and cold fluid were smoothed to 
match the laminar flame thickness. 



We consider two types of flows, sinusoidal shear flow 
with amplitude U and wavelength L, perpendicular to 
the initial temperature front (s), 



u = U I 0, cos 



2lTX 



(2) 



and cellular flow with amplitude U and wavelength L, 

/ . 2irx 2ny 2ttx . 2iry\ 

u = u I sin — — cos — — , — cos —j— sin — — I . (3) 
V Li hi Li Li J 

In the quenching simulations, the size of the cell L/2 was 
a fraction of W, so that the initial band always contains 
integer number of cells. 

Most of the reaction enhancement computations were 
done using KPP reaction rate [^7| |3(| in the advection- 
reaction-diffusion equation (|l|), 



R(T) = ~T(1-T), 



(4) 



with some of the simulations repeated with ignition type 
reaction, 



R(T) 



T 



(1-To 



(i-n 



(5) 



where T represents threshold temperature, below which 
R(T) = 0. In quenching studies we use ignition type re- 
action (||) with threshold temperature T — 0.5. Both 
reaction rates (||) and (||) were chosen to exactly match, 
in the absence of advection, laminar burning velocity 
i>o = \J k/t. The corresponding laminar flame thick- 
ness is in both cases of the order of 5 = y/nr; for KPP 
reaction it is several times wider than for ignition. 

Equation (|l|) with reaction rates (f|) and (||) has 
been solved using a fourth-order explicit finite difference 
scheme in space and a third-order Adams-Bashforth in- 
tegration in time. The grid size was chosen so as to ac- 
curately represent the shear across the reacting region: 
typically of the order of 12 zones across the flame inter- 
face for thin fronts and at least 32 zones per period for 
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FIG. 1: Schematic representation of initial conditions and 
velocity field. Problem setup, from left to right: reaction 
enhancement in shear flow; reaction enhancement in cellular 
flow; quenching in shear flow; quenching in cellular flow. Dark 
tone corresponds to the cold fluid (T = 0) , light - to the hot 
fluid (T= 1). 
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thick fronts. The computational domain extended a con- 
siderable distance upstream and downstream from the 
burning front so that boundary effects were negligible. 
In flame enhancement simulations the overall grid was 
remapped following the propagation of the front, thereby 
allowing for long integration periods — of the order of 
1000 reaction times r. We found that these long integra- 
tions were necessary in order to reproduce correctly the 
asymptotic behavior of the propagation speed in the case 
of strong advection. 

As a measure of the reaction enhancement we use the 
bulk burning rate 



v(t) = 







dT(x,y,t) 



at 

R(T) dydx 



dydx 



(6) 



The second equality in (g) can be justified by integra- 
tion by parts. The bulk burning rate coincides with the 
front velocity in a case where the solution is a travel- 
ling wave, but provides a more flexible measure of com- 
bustion. Physically, v(t) can be understood as the total 
amount of reacted material or the total heat production. 
In all simulations done for reaction enhancement, v(t) 
approaches an asymptotic value (for cellular flows one 
should average in time to arrive at this value) and we 
denote this asymptotic value v. 

We remark that in shear and cellular flows equation 
(|l|) admits travelling wave-type solutions, called pulsat- 
ing fronts (see e.g. & 10, 0,0, ra). A rigorous stability 
theory for these solutions exists but is not complete (see 
p3| for a recent review). Except when quenching occurs, 
in our simulations we always observed convergence of the 
solution to such waves, so the bounds on v also provide 
bounds for the propagation speeds of pulsating fronts. 

In the quenching studies, we measure the total amount 
of burned material per wavelength, 



T(x,y) dy dx, 



(7) 



which can also be interpreted as the width of the non- 
perturbed horizontal band with temperature T = 1 sur- 
rounded by the fluid with T = Q. This quantity is related 
to the bulk burning rate per interface, v(t) = w(t). De- 
pending on initial conditions and flow parameters, w(t) 
either approaches an asymptotic value, e.g. v(t) — > 0, or 
increases with constant rate, that is v(t) — ► v (for cellu- 
lar flow, in the time-averaged sense). In the first case we 
say that flame quenches; the main objective of quenching 
simulations is to determine under which conditions this 
happens. 



III. SHEAR FLOW: REACTION 
ENHANCEMENT 

We carried out simulations with the sinusoidal shear 
flow with amplitude flow U , and wavelength, L, given by 
Eq. (|J). In this section we are interested only in reac- 
tion enhancement phenomena, and therefore for initial 
conditions we consider T = 1 in the semi-infinite domain 
y < 0, and T = for y > 0. We carried out computa- 
tions for both KPP and ignition type reactions, but did 
not find significant differences in the qualitative behav- 
ior. The numerical results presented in this section are 
obtained with KPP reaction term (^). 

Of special interest is the dependence of the effective 
propagation rate v on the velocity amplitude, U, and 
wavelength, L, which defines the characteristic length 
scale of the flow (Fig. ||). For small amplitudes, U <C v , 
our results are in agreement with the quadratic law 
v ~ v + elf 2 , which goes back to Clavin and Williams 
Jl7[ for turbulent flow and has been recently proved rig- 
orously for shear flows in p9fl . We did not study this 
regime in detail since our main interest is in the strong 
advection case. 

For the amplitudes U 3> v the results are in good 
agreement with the linear law v = all + b, where the 
coefficients a and b < v depend on the geometry of the 
flow. In the situation where the scale of the flow is much 
larger than the reaction length scale, L 3> S, our data 
agrees with v = U + v a ■ This law has been proposed in 
Q for shear flows which vary slowly compared to the 
typical reaction length and rigorously proved in [ pp[ un- 
der similar assumptions. For any fluid flow, the regime 
L ^S> S is closely related to the so-called geometrical op- 
tics combustion regime E3] , the limit where reaction time 




FIG. 2: Bulk burning rate ([]) as a function of the shear flow 
amplitude for different shear wavelengths. 
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FIG. 3: Isotherms within the front in the geometrical optics 
limit. Here L/8 = 256, and U/v D = 4. 



and length scales approach zero. In the framework of the 
equation (^) this corresponds to the limit k, t — > while 
k/t remains constant. Indeed, by rescaling equation (^) 
with a factor L/8 in space and time, we find that the bulk 
burning rate v for the original equation ([!]) is the same 
as for equation with modified diffusion and reaction (pj) 



S L 

T t + u-VT- k-AT = — R(T). (8) 

L OT 

As L grows, equation (Q) approaches the geometrical op- 
tics limit. 

Quite often, front propagation in the thin front and 
fast reaction limit is modelled by Hamilton-Jacobi type 
equations. One such model is the G-equation 



G t + u-VG = v \VG\, 



(9) 



where the front is defined by a constant level surface of 
the scalar G (see e.g. Q). The G-equation describes 
propagation of the front according to Huygens principle; 
that is, the front (i) is transported by fluid flow, and (ii) 
propagates normal to itself with the speed v Q ■ The law 
v = U + v can also be understood from the point of 
view of geometrical optics since it is easily derived from 
the G-equation. Recently Majda and Souganidis pointed 
out that the G-equation does not provide the geometri- 
cal optics limit of the reaction diffusion equation (^) in 
a precise sense |25|, [$7|. However the rigorous bounds 
derived for the true effective equation still give the same 
prediction for v in the case of shear flow [p7| . 

In the situation where L becomes comparable to 5, 
the coefficient proportionality a between velocity ampli- 
tude and flame propagation rate is no longer equal to 
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FIG. 4: Bulk burning rate in high wavenumber sinusoidal 
shear flow as function of shear amplitude (points), compared 
with isotherm elongation (solid lines) given by Eq. 



unity. The rigorous lower bound for v from [p5[ takes 
form v > G\U 1+ g , where n — 2ir5/L. This bound is in 
good qualitative agreement with an argument proposed 
by Abel, Celani, Verni and Vulpiani jjj based on the ef- 
fective diffusivity for the shear flow. It is well known that, 
if the problem is considered on sufficiently large time and 
length scales, the effect of the advection of passive diffu- 
sive scalar can often be modelled by effective diffusivity 
0, H . The expression of effective diffusivity in a strong 

shear flow goes back to Taylor fjsfl , Keg = K +5 f ^Tm) K ' 

where m = 2ttS/1 and / is the typical length scale of the 
flow. In the presence of reaction, we take I = min(<5, L), 
since the advection balances with reaction instead of dif- 
fusion if L > (5. This leads to the qualitative prediction 
v ~ U if L > S and v — UL/S if L < S. We obtained 
good although not perfect agreement with this predic- 
tion. This is not surprising given the heuristic derivation 
of the expression for the effective diffusivity and its pos- 
sible dependence on more subtle geometric properties of 
the flow. 

Additional understanding of linear dependence v(U) 
can be gained from studying the relationship between the 
burning enhancement and the structure of the front, in 
particular level sets of the temperature. Assume that in 
the geometrical optics approximation the front is given 
by the function y = /(x); then for the travelling wave 
obeying Huygens principle and propagating with speed 
v, we have 



where u(x) is profile of the shear, u[x) — Ucos^-. In 
the case where u is a mean zero flow, this leads to the 
expression 



(10) 



Thus, we obtain a well-known fact that the speed of prop- 
agation is proportional to the area of the front which in 
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L/8 = 32 

T = 0.001 




FIG. 5: Isotherms within the front for cases with L/5 = 
32 and U/v = 10,12,14. The individual curves have been 
rescaled by a factor of (C//i> ) -1 in the y direction. 

geometrical optics limit coincides with a level set of T 
(see Fig. || for a picture of level sets in a situation close 
to geometrical optics). It is interesting to test to what 
extent this relationship remains true in the situations 
where geometrical optics regime is no longer valid, for 
example when L is comparable to S. We found that for 
large U, there is still good agreement between the elon- 
gation factor of the level sets of temperature and com- 
bustion enhancement (Fig. ^). Moreover, we found that 
for large U the temperature distribution across the front 
scales with the shear amplitude (see Fig. |^), providing 
another explanation of the linear dependence of the bulk 
burning rate on the amplitude of the flow. This scaling 
behavior can be understood in terms of the approximate 
self-similarity of equation with respect to the change 

of variables y — jjt^;, x = x/S, t = t/r, which gives 

The only term which now depends on U is the one pro- 
portional to the second derivative in y, and it becomes 
negligible in the limit of large U. Indeed, equation ( |ll| ) 
without this term is hypoelliptic, and so addition of the 
second derivative term does not constitute a singular per- 
turbation. That leads to [/-independent propagation rate 
v = jj — v(L) and to linear proportionality v oc U for 
large U/v . 

We conclude this section by a remark that understand- 
ing of the combustion enhancement in a shear flow ap- 
pears to be useful in some situations where flows with 
different structure are involved. In particular, in the 



reactive Boussinesq system the flow consists of vortices 
moving along with the reaction front [j54j| . However in 
the frame moving with the front the effect of such vor- 
tical flow is similar to the shear. The prediction for the 
reaction enhancement in such system based on the re- 
sults for the shear flows appears to be in agreement with 
numerically observed behavior pm . 



IV. CELLULAR FLOW: REACTION 
ENHANCEMENT 

Cellular flows have been studied by many authors (see 
e.g. p"5| , |22j), since similar fluid motions appear in many 
important applications; classical examples are the two- 
dimensional rolls of the Rayleigh-Benard problem and 
Taylor vortices in Couette flow. 

For cellular flow simulations, we use a velocity field 
given by Eq. (|^). The flow is controlled by two param- 
eters, velocity amplitude U and wavelength L. As in 
simulations of reaction enhancement by shear flow, we 
consider initial conditions with T = in the upper half 
of the computational domain (y < 0), and T — 1 in the 
lower half (y > 0). Most of the results presented in this 
section were obtained with KPP reaction term (Q); the 
influence of the reaction type on the reaction enhance- 
ment is discussed at the end of this section. 

There are several regimes, which can be classified ac- 
cording to the relations between the characteristic scales 
present in the problem. There are three characteristic 
time scales: advection, tjj = L/U; reaction, tr — t; and 
diffusion, tjj = L 2 /k. In this paper, we mostly stud- 
ied two regimes: tjj <C tr <C Td and tr <C tjj -C tjj. 
The first is the regime of strong advection; in the sec- 
ond regime advection can be very strong as well, but 
is compensated by large cell size L, so that the reac- 
tive time scale becomes the fastest in the problem. If 
Td <C tr, or equivalently, L <§; <5, we have diffusive or 
small cell regime. The remaining regime corresponds to 
tr <C td "C tjj, and therefore [/ « d . Hence we have 
slow advection; this situation is of less interest to us since 
v < U + v under very general conditions p9| , so the ef- 
fect of the advection on the propagation speed is minor. 

In the regime tr <C tjj <C tjj our simulations show 
good agreement with geometrical optics models. In Fig- 
ure we show a typical picture of flame in the regime 
close to geometrical optics. The G-equation, (^), which 
is closely related to geometrical optics regime, is invari- 
ant under simultaneous rescaling of time and space by 
the same factor. This suggests that for flames approach- 
ing the geometrical optics limit, one should observe this 
similarity, and indeed we do as shown in Fig. ^. The 
earliest prediction of the front propagation speed in a 
cellular flow for large U within the geometrical optics 
framework appears to be due to Shy, Ronney, Buckley 
and Yakhot |47j. Using heuristic reasoning, they pro- 
posed that v ~ U/ \og(U/v ) if one considers front ad- 
vancing according to Huygens principle. The same law 
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FIG. 6: Burning regimes in cellular flow: (I) tr < to < 
tu — slow advection geometrical optics; (II) tr < tu < td — 
fast advection geometrical optics; (III) tu < tr < td - 
fast advection with radial burning within cells; (IV) tu < 
td < tr — fast advection with uniform burning within 
cells; (V) td < tu < tr — limited advection in small cells; 
(VI) td < tr < tu — slow advection in small cells. 



is proposed in [9 based on more detailed analysis. To 
illustrate the origin of this law, let us sketch an argu- 
ment providing the lower bound for v. Let us look at the 
propagation of the flame tip along the path ABCDE on 
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q. Assuming that at every point of ABCDE the 
ocity is given by u{x, y) + v , and integrating in 
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FIG. 7: Flame in cellular flow with amplitude U/v a = 20 
and period L/8 = 1024 (upper row) and L/S = 512 (lower 
row) . Snapshots for the first case were taken with time inter- 
val 24 8/vo and for the second - with time interval 12S/v - 
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This lower bound, when doubled, gives a very good fit 
to our numerical data and is represented by a solid line 
m Fig. | As one can see in Fig. 0, the tip of the flame 
follows the path close to ABCDE, but avoiding corners; 
this may account as a factor for the difference in speed 
compared with the lower bound. Our results in the ge- 
ometrical optics regime agree with results of ||; how- 
ever we stress that our simulation was carried out for the 
reaction-diffusion equation (nh using the same numerical 
scheme in all regimes, while [fa] uses a different model in 
geometrical optics limit. We remark that in [jl], ^] the 
possibility of the regime, where tr <C to and v behaves 
as t/ 3 / 4 , was proposed. We could not definitively con- 
firm existence of such a regime due to the closeness of 
(U/v ) 3 / 4 and (U / v ) / log({// v ) curves in the range of 
tested parameters (Fig. ||). 

It should be emphasized that the geometrical optics 
regime requires not only the thin front assumption, L 3> 
5, but also fast reaction in comparison with advection, 
tr <C tu; in other words, velocity must be limited by 
U <C (L/5)v Q . When this restriction is broken we observe 



significant decrease in flame propagation speed compared 
to the geometrical optics prediction (Fig. [)]). Figure [h] 
further illustrates this point, showing that on a logarith- 
mic scale, there is a marked change in the slope of {v/v ) 
as a function of U as tr/tjj increases. When t\j exceeds 
tr, we observe power-law, v ~ U 1 ^ 4 , as proposed by Au- 
doly, Berestycki and Pomeau Q , and confirmed in Q in a 
narrower range of parameters. The measurement of the 
slope 1/4 was sufficiently precise to distinguish it from 
the v ~ J7 1 / 5 behavior, a lower bound rigorously proved 
in |3^|. The observed v ~ [7 1 / 4 scaling extends to the 
limit of cell sizes small in comparison with the laminar 
thickness, L <C S. We remark that the laminar front 
thickness for the KPP reaction is of the order of 16 5 
which is large compared to smallest cell size L/2 = 4 6 
shown in Fig. QG; limited data available for L/S = 4,2 
(shown in Fig. |ll|) also confirms the v ~ U 1 ^ 4 scaling. 
In the very small cell regime, the v ~ J7 1 / 4 scaling was 
rigorously proven in ^| using homogenization approach. 

We also studied the dependence of v on the cell size 
while XJ / v is fixed. Figure |l2| shows changes in the struc- 
ture of the flame with the increase of the cell size — from 
a more diffusive front to a front approaching geometrical 
optics behavior. The flame propagation speed, normal- 
ized by (U/vo) 1 / 4 factor, is presented in Fig. jn|, As L/S 
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FIG. 



Approximate path of the tip of the flame 



increases, we see the transition from the tjj -C tr regime, 
to the geometrical optics, where flame propagation speed 
is independent of cell size. For tjj <C tr the data collapse 
to a single curve, suggesting the power scaling with L/S, 
with power changing from 1/4 for small L/S to 3/4 for 
large L/S. The resulting scaling can be summarized as, 



v/v 

V/V 

v/v 



(U/v ) I \0g(U/v ), T R « Tjj « T D , (13) 

(U/v ) 4 / 4 (L/S) 3 / 4 , r [/ «r ii «r D . (14) 
(U/v ) 4 / 4 (L/S) 1 / 4 , r [/ «r D «r fl , (15) 



To explain observed the flame propagation speed, let 
us consider a model based on the effective diffusivity, 
proposed by Audoly, Berestycki and Pomeau ||. When 
velocity is high, tjj <C tr, the sharp temperature gra- 
dients appear in the narrow boundary layer at cell bor- 
ders (Fig. |l3|). The thickness of the boundary layer, h, 
is determined by the balance of diffusion and advection, 
n/h 2 — U / L, and is much smaller than 5, 



y/nL/U. 



(16) 
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FIG. 10: Flame propagation velocity as function of the ratio 
of laminar burning time to the vortex turnover time. The 
L/S = 2048 prediction is based on doubled Eq. (^) for geo- 
metrical optics. 



(This argument is a concise, naive version of the con- 
siderations appearing in the derivation of the effective 
diffusivity in cellular flow, @, [Sj, |j| f||.) 

A discrete diffusion equation modelling the original 
equation (0) has been suggested in @, ||] : 



at L l t 

Here 9 n is the average temperature in the n-th cell; and 
K e ff is effective diffusivity 



K c ft = nL/h. 
This leads to the propagation rate, 



(18) 



v = \J k c «/t ~ v a y/L/h, 

and, taking into account Eq. (|l6|), to the scalin g ( 15[ ) . We 
found that the prediction for speed given by (J17 ^agrees 
with numerical simulations only if L < 5. 




FIG. 9: Flame propagation velocity in thin flame regime as 
function of the cellular flow amplitude. Solid line is doubled 
Eq.@, dashed line is v/v a = {U/v ) 3/4 . 
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FIG. 11: Dependence of the flame propagation velocity on 
the size of the vortex. 
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FIG. 12: Flame in a cellular flow with amplitude U/v = 100 
and period L/S = 16, 64, 256 (left to right) 



For large cell sizes, L^> S, the diffusive model (17) no 
longer accounts fully for the flame propagation process. 
The main objection to the model is that (|l7|) assumes 
temperature is uniform inside the cell, and the reaction 
term can be estimated at the average cell temperature. 
However the numerically observed behavior demonstrates 
at first sharp temperature gradients at the border of the 
cell, later evolving into flame propagation inside of the 
cell roughly at the laminar flame speed (see Fig. [D] for 
the structure of the front inside cells). Indeed, the cel- 
lular flow has no efficient mechanism for mixing between 
the streamlines, and the diffusion time scale in that di- 
rection is of the order L 2 /k (practically not enhanced) 
EJ. Therefore the combustion process inside the cell 
takes time of the order of L/v Q , rather then S/v (which 
corresponds to substituting the average temperature into 
the reaction term) . 

Here, we will modify the effective diffusivity model (|i~7|) 



to account for slower burning inside large cells. As in the 
case of small cells, the flame propagation from one cell 
to another is enhanced because of the high temperature 
gradient in the boundary layer with width given by ex- 
pression (16|). The heat coming to the cell through the 
cell boundary is distributed on the scale of 8 (as opposed 
to L, in the case of small L). We further notice that 
due to the fast advection, the temperature is essentially 
equal along the streamlines inside the cell (see Fig. [l2| ), 
which allows the flame to propagate directly from one 
boundary layer to another (Fig. [l3|). That allows us to 
write the discrete diffusion equation similar to (|j~7j), but 
replacing averaging in the cell by the averaging in the 
strip of width S along the cell border ((S-layer), 



dt 



Kcff 

s 2 



[#71-1 



29 n + 9„ 



-R(0 n 



(19) 



Here 6 n is the average of the temperature in the (S-layer 
of the n-th cell; and K c g is effective diffusivity, 



K c g = nS/h. 



(20) 



The temperature in a laminar front varies on the scale 
5, and so estimating reaction in a <5-layer by R(9 n ) is 
justified. Equation (|l^) does not take into account the 
heat flux from the i5-layer to the bulk of the cell. How- 
ever, since h <§; S, this heat flux does not enter the 
main balance. In other words, the front described by 
( ]l9| ) propagates entirely in the 5-layers, with the rate 
vs = \Zk c s/t ~ v ^jS/h. Substituting h from Eq. ( |i"6| ) 
we obtain, 



v s ~ vo (U/voY^iL/S) 



-1/4 



(21) 



The time needed to ignite a new cell, e.g. to warm up 
a layer of the size of the order S in that cell, is equal to 
ts = 5/vg. Once the width of a warmed up layer reaches 
size of the order of S, reaction becomes capable of sustain- 
ing the temperature. Further propagation of the flame 
corresponds to the basically laminar front movement in- 
side the cell, and takes time r CG n ~ L/v Q . 

The total bulk burning rate v is of the order of v times 
the number of burning cells, which can be estimated as 
the ratio of the cell burning time, r ce n ~ L/v , to the 
time needed to ignite a cell, ts = S/vs- Therefore the 
bulk burning rate is obtained by multiplying the number 
of burning cells with v , or essentially by normalizing vg 
with a factor of L/S, 



Tccll L 

V — V ~ — VS ~ Vo {U/Vo 

TS 



,V4 



(L/S) 



3/4 



which agrees with numerically observed scaling (14). 

Finally, we would like to mention the effect of the re- 
action rate. Similar to the shear flow, we found that the 
reaction type does not influence asymptotic scaling laws 
like (|l^) or (|l4|), although there is certainly a difference 
in the constant factors. However, in cellular flows there is 
an interesting phenomenon which is present for ignition- 
type but not KPP reactions. The dependence v(U) is 
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FIG. 13: Temperature profile at the middle of the cell (solid) and temperature averaged in x-direction (dashed) for L/S = 
8, U/vo = 64 (left) and for L/S = 256, U/v a = 128 (right). Bottom plots are blow up versions of top plots. Dotted line in the 
top plots represents laminar front stretched by factor v/v a in {/-direction, while in the bottom plots it represents non-modified 
laminar front. 



always monotone increasing in the KPP case, but it may 
exhibit a temporary reversal for ignition-type reactions 
(Fig. |l4|). This effect has been discovere d by Kagan and 
Sivashinsky ||(| and further studied in |[3lf . We found 
that this phenomenon is more pronounced when the re- 
action threshold T a is closer to unity, in agreement with 
the arguments of plf. 



KPP 

ignition T = 3/4 
ignition T = 1/2 
ignition T = 1/4 




1 2 4 



16 32 64 128 256 512 
U/v„ 



FIG. 14: Flame propagation velocity for different reactions 
(L/S = 8). 



V. QUENCHING 

In this section we address another effect that advec- 
tion can have on the combustion process — quenching. 
We say that reaction is quenched if the average tempera- 
ture goes to zero uniformly with time. Quenching occurs 
in the systems with ignition type reaction when, due to 
diffusion and advection, temperature drops below the ig- 
nition threshold everywhere and the integrated reaction 
rate becomes identically zero. 

If the size of the region is small enough quenching 
can be caused by diffusion alone, e.g. without advec- 
tion, as shown by Kanel |32|. Kancl considered the 
one-dimensional reaction-diffusion equation T t — nT xx — 
r _1 i?(T). He found that there exist two critical sizes 
Wo < W* such that if the initial size of the hot region 
(where T = 1) is smaller than W a , reaction quenches, 
while if the initial size of the hot region is greater than 
W* , two fronts form and propagate in opposite direc- 
tions. Two different critical sizes are likely an arti- 
fact of the proof; in our simulations, we always found 
W = W* ~ S and will refer to the single critical size W . 
In two and three dimensions, the presence of advection 
may lead to stretching of the initial hot spot, thus mak- 
ing diffusion more efficient at cooling, and consequently, 
at quenching. 

In our numerical simulations we study quenching un- 
der the influence of advection, in particular in shear and 
cellular flows. As in previous sections, the prescribed 
flow velocities are defined by Eq. @ for shear flow and 
by Eq. (||) for cellular flow. For all simulations we used 
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FIG. 15: Sequence of snapshot of temperature distribution in the shear flow with L/8 — 4 (top) and in the cellular flow 
with L/8 = 4 (bottom). Initial condition was a hot band of the width W/8 = 6 and W/5 — 4 for shear and cellular flow 
simulations respectively. Velocities amplitudes are below critical on the left and above critical on the right (shear: U/v ~ 13 
and U/vo = 14; cellular: U/v a = 600 and U/v = 800). The time is given in units of r. 
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FIG. 16: Temperature integrated over area < x < L, 
— oo < y < co per wavelength L for different values of shear 
amplitude U, measured for L/8 — 6 and W/S — 6. 



the ignition- type reaction (||) with threshold T Q = 1/2, 
since quenching cannot occur for the KPP-type source 
term [[l3[ Q. As the initial conditions we use a hori- 
zontal band of width W with temperature above critical 
(T = 1) within the band and with temperature below 
critical (T = 0) outside the band. 

The typical evolution of the system described above is 
shown in Fig. Il5[ We found that the temperature distri- 
bution in both for shear and cellular flows evolves accord- 
ing to one of the two possible scenarios, depending on the 
amplitude of the advection velocity. For lower advection 
velocities, after an initial transient period, the system de- 
velops solution characterized by a wide, steadily growing 
burned region between two wrinkled fronts propagating 
in opposite directions. These fronts are exactly as de- 
scribed in the preceding sections with regard to structure, 
speed, and dependence on the flow properties U and L. 
For higher advection velocities, the temperature eventu- 
ally drops below T everywhere, after which no burning 
occurs. We denote by U CI the value of advection veloc- 
ity which triggers the system between these two scenarios 
(further we refer to them as burning and quenching) , and 
establish the relation between U CI and the initial condi- 
tions and structure of the flow. 

We measure the critical value of velocity amplitude 
U CI (L,W) using the following procedure. For each com- 
bination of the initial hot band size W and velocity pe- 
riod L, we execute a number of simulations for different 
velocity amplitudes U. For each simulation we measured 
the total amount of burned material per period w(t), de- 
fined by Eq. (@), as a function of time (shown in Fig. |16|). 
The burning systems (with higher velocities, where two 
front are formed) are characterized in Fig. [l^ by constant, 
non-zero slopes, corresponding to constant reaction rate. 
These rates are independent of initial conditions, and 
are equal to double the burning rate v(U, L), since there 
are two fronts. The quenched systems are characterized 
by evolving to constant w(t). Both formation of steady 



FIG. 17: The critical amplitude of shear flows with different 
wavelengths as function of initial width of hot band. 



fronts and quenched solutions requires some transition 
time. 

The summary of results for a sinusoidal shear flow is 
given in Fig. [17| and Fig. Jig . We found that £/„■ scales 
linearly with W (see Fig. p7^, 



U r , = a- 



W 



(22) 



as predicted in j20|, with the coefficient a strongly de- 
pendent on the wavelength of the advection velocity 
(Fig. |l^). Shear flow is most effective at quenching in 
the intermediate range of wavelengths, namely when L is 
of the order of a few reaction lengths 5. The quenching 
mechanism for small and for large L is different; one has 
to distinguish between the ability of the flow to stretch 
the front over the larger scales and to make the initial 
hot band uniformly thin. 

For small L, the rapid spatial variation of the flow 
velocity is well approximated by effective diffusion. The 
effective diffusivity for strong shear flow scales as k c s ~ 
U 2 1 2 /k, where / = mm(S, L). The characteristic length 
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L/5 



FIG. 18: Dependence of the factor a in Eq. (|22j) on the wave- 
length of shear flow L. Measurements were taken at W/S — 6 
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FIG. 19: Growth in the value of [/ cr for a fixed W as the size 
of fiat plateau increases. 

scale for reaction with this renormalized diffusion behaves 
as l e ff ~ y/K c gr and scales linearly with U. Then Z e ff ~ W 
gives U CI ~ jWt~ 1 for small L. We remark that in the 
limit of small L, the effective diffusivity argument can be 
justified by a rigorous homogenization procedure p9| . 

For large L, the nature of quenching is related to the 
appearance of the (almost) constant regions in the veloc- 
ity profile. We observe the behavior U CI ~ L 2 for large 
L, which can be explained in the following way. In order 
for quenching to occur, the shear flow should stretch the 
initial hot region thinner than Kanel's critical length (of 
the order 6) in time less than the reaction time r, so that 
reaction does not have time to compensate cooling by 
advection. The stretching is least efficient near the tip 
of the velocity profile. At the tip of a sinusoidal profile, 
the difference between flow velocities at two points sepa- 
rated by a distance 8 is U(5/L) 2 . Therefore we obtain a 
sufficient condition for quenching U cr t ~ W, which 
leads to 

We also examined a degenerate case of shear flows with 
a plateau in the velocity profile. For such flows it has 
been shown in |p0[ that quenching does not happen as 
soon as the size of plateau is larger than certain critical 
size of the order 8 and the size of initial bandwidth W 
exceeds W ■ As expected, we found that U cr diverged to 
infinity as the size of the plateau approached a critical 
value (Fig. [n]), in agreement with results of [^o| . This 
phenomenon can be understood in terms of reaction and 
diffusion alone: in the region where the profile of the ve- 
locity is flat, there is no stretching of the initial hot band. 
If the size of the hot band is roughly larger than Kanel's 
critical size W , then reaction can compete with diffusion, 
there will be no quenching, and eventually propagating 
fronts will form. 

Quenching in cellular flow requires significantly higher 
advection amplitudes. For relatively small cell sizes 
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FIG. 20: The value of U c t for cellular flows with different 
periods. 



L < 8 where quenching is possible, we find that the crit- 
ical velocity U cr satisfies U cr ~ W 4 (see Figure ^). No- 
tice that this correlates with the dependence v ~ U 1 ^ 4 for 
the speed of advection-enhanced flame. This correlation 
is not coincidental, and can be explained as follows. The 
speed up law indicates that the size of the region where 
reaction happens in a stabilized combustion regime scales 
as {U/voY^iL/Sf^S for large U. If the width of the ini- 
tial band of hot material is of the order smaller than 
(U /voY^iL/Sy^S, advection carries away the energy of 
the hot material over the larger region faster than reac- 
tion is able to compensate the falling temperature. This 
leads to quenching. 

We remark that quenching is impossible if the cell size 
is sufficiently large, L > 8, and W > 8. The reason is 
similar to the flat plateau effect in the shear flow. Fluid 
advection does not provide mixing inside cells in the di- 
rection perpendicular to the streamlines, and thus if the 
cell is large enough reaction can sustain itself against dif- 
fusion. This result has been proved in Jy] . 



VI. CONCLUSIONS 

We carried DNS calculations of an advected scalar 
which reacts according to a nonlinear reaction law. We 
studied combustion enhancement and quenching phe- 
nomena in two typical classes of flows, shear and cellular. 
In a shear flow, we find linear dependence v — all + b of 
the combustion speed v on the amplitude of the flow U 
in the strong flow regime. The factor a depends on the 
relationship between the period of the flow L and typical 
reaction length scale 8, is equal to 1 if L ^> 8 and tends to 
zero if L/8 — * 0. The observed behavior is in agreement 
with recent rigorous |l9[ and numerical |j| results. In a 
cellular flow we studied primarily two regimes character- 
ized by the relationships td < W < tr and Tp < tr < 
tjj between diffusion, reaction and advection time scales. 
We found that combustion speed in the first regime is 
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close to predictions of models based on geometrical op- 
tics limit, v ~ XJ I \og(U/v ). In the second regime where 
large U dominates, we found v ~ v (U /v ) 1 / i (L/5] 3 / 4L . 
This agrees with the prediction of the effective diffusion 
model 0, ||, [| in terms of the power of U but has different 
dependence the cell size L. We proposed an explanation 
of the observed behavior with a modified effective diffu- 
sion model where enhanced diffusivity is concentrated in 
the boundary layers. 

As opposed to combustion enhancement, quenching 
may happen if the reaction term is of ignition type and 
initial temperature is higher than critical in a finite re- 
gion. If the shear flow velocity profile does not have a 
plateau of sufficiently large size, or the size of the cells in 
cellular flow is not too large, then for any initial hot band 
size W there exists U CI (W) such that for U > U cl - quench- 
ing takes place. If U < U CT , two fronts form and propa- 
gate with the speed of the developed advection-enhanced 



front. In the case of shear flow, U c[ depends linearly on 
W with a factor a(L)/r. Quenching is most efficient for 
the flows with L on the order of a few typical reaction 
lengths 6. For the cellular flows, U cr scales as W 4 . The 
results are in good agreement with theoretical arguments. 
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